Entanglement growth and simulation efficiency in one-dimensional quantum lattice systems 
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We study the evolution of one-dimensional quantum lattice systems when the ground state is perturbed by 
altering one site in the middle of the chain. For a large class of models, we observe a similar pattern of entan- 
glement growth during the evolution, characterized by a moderate increase of significant Schmidt coefficients 
in all relevant bipartite decompositions of the state. As a result, the evolution can be accurately described by a 
matrix product state and efficiently simulated using the time-evolving block decimation algorithm. 

PACS numbers: 



I. INTRODUCTION 

The numerical study of quantum many-body systems is a 
challenging computational task due to the exponential growth 
of the Hilbert space dimension with the system's size. Lat- 
tice systems in one spatial dimension, such as quantum spin 
chains, are a noticeable exception. There, the density ma- 
trix renormalization group (DMRG) algorithm allows for the 
precise computation of ground state properties |0], while 
the time-evolving block decimation (TEBD) algorithm can be 
used to simulate time evolution [00]. Such techniques pro- 
vide valuable insight into quantum systems, therefore facil- 
itating progress in several forefront areas of research, both 
in science — e.g. condensed matter, quantum optics, atomic 
and nuclear physics, quantum chemistry — and technology — 
e.g. quantum information processing, quantum computation, 
nano-technology. 

Many-body entanglement is at the very core of the achieve- 
ments of the DMRG and TEBD algorithms. In both cases, 
a key ingredient is the use of a matrix product state (MPS) 
to represent the state \^) of the system fl]. A MPS leads to 
an efficient representation of |^) provided that the amount of 
entanglement in the system is sufficiently small. Thus, the 
success of the DMRG relies on the remarkable fact that in 
one spatial dimension the ground state of most local Hamilto- 
nians has only a limited amount of entanglement. Likewise, 
the TEBD algorithm, based on updating the MPS in time, is 
efficient as long as no large amounts of entanglement are pro- 
duced during the simulated evolution. 

Identifying time evolutions that only involve small amounts 
of entanglement is, consequently, of central importance in or- 
der to establish the range of applicability of the TEBD al- 
gorithm — and of a whole score of subsequent proposals 
fi 0.0] > including implementations within the DMRG for- 
malism la 0,111] often referred to as tDMRG, that are based on 
the same idea, namely on adapting the MPS representation of 
the state of the system, so as to account for its changes during 
a time evolution. At the same time, a better understanding of 
the dynamics of entanglement in one-dimensional many -body 
systems [9] is of interest in the areas of quantum informa- 
tion processing, where entanglement is regarded as a useful 
resource, and of quantum computation, where entanglement 
is necessary in order to achieve a significant speed-up with 
respect to classical computers Bfl^l. 



In this paper we study the generation of entanglement in 
a particularly relevant class of time evolutions. Specifically, 
we consider a chain of TV sites, where each site is repre- 
sented by a finite-dimensional Hilbert space and labeled by s, 
(s = 1, • • • , TV). The system, initially prepared in the ground 
state I^gs) °f some local Hamiltonian H, is perturbed at time 
t = by applying an operator A on lattice site sq. As a re- 
sult, the state of the system 1^(0)) = A s ° \^ as ), no longer an 
eigenstate of H, evolves non-trivially in time, see e.g. Fig. 
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During this evolution entanglement is produced, adding to 
the entanglement that might already be present in the ground 
state. 

Notice that state (^(t)) appears in the computation of the 
unequal time correlation function 



(B s {t)A s °(Q)) = (y GS \e lHt B s e~ lHt A B0 \ty o 
= e lEot (y os \B s \y(t)), 
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where a second operator B has been applied on site s at time 
t. Therefore, by characterizing the growth of entanglement 
in evolution (Q~|i, we will be able to assess the efficiency with 
which the algorithms of Refs. [2, 3} 0] can be used to 
compute the two-point correlator ©. We recall that from this 
correlator one can extract quantities such as Green's functions 
or dynamic structure factors, and thus learn about a number 
of properties of the system, including its response to external 
probes, e.g. neutron or photon scattering. 

The study of the entanglement and simulatability of a time 
evolution of the form ([T|i was initiated in Ref. I3J and is in- 
timatedly related to the development of the TEBD algorithm. 
Originally, the TEBD algorithm was created in order to char- 
acterize the role of entanglement in quantum computation 10] ■ 
Soon afterward it was noticed that, for systems in one spa- 
tial dimension, low energy time evolutions such as (Q~|i seemed 
to often involve small amounts of entanglement [0], precisely 
in the way that would allow the TEBD algorithm to work effi- 
ciently. This opened up the possibility to simulate this class of 
dynamics. One may claim, a posteriori, that it is somehow ex- 
pected that no much entanglement will be created when a sin- 
gle local perturbation is introduced on the ground state \^ as } 
of H. However, the subject is far from being well understood 
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and, for instance, given H and local operator A, there are no 
known conditions to guarantee that the evolution can be sim- 
ulated efficiently for reasonably long times. The best result in 
this direction, due to Osborne fTTIl , shows that that time evo- 
lutions in one spatial dimension can be simulated efficiently 
(that is, with polynomial resources in N) only up to small 
times t w log(iV). Instead, numerics suggest that the simula- 
tion of ([T]i can be performed efficiently for times that scale at 
least as the size of the system, t ~ N, JHHt]. 

Several authors (see for instance [HI [II Bill]) have 
already considered particular instances of evolution ensuing 
either a local perturbation of the ground state as in Eq. (Q~|), 
or a local quench (which is qualitatively equivalent to a lo- 
cal perturbation). In these specific cases it has been observed 
that, indeed, the amount of entanglement remains sufficiently 
small as to allow practical simulations for large times and sys- 
tem sizes (in a sense further specified in the next section). The 
goal of the present paper is to establish how general this re- 
sult is, by conducting a systematic numerical study of evolu- 
tions of the form (HJ in one-dimensional models, including: (i) 
systems of spins, fermions and bosons; (ii) gapped and gap- 
less phases; (in) homogeneous and disordered interactions; 
and (iv) integrable and non-integrable systems. We find that, 
whereas for each model the evolution may represent a very 
different physical process, in all cases the dynamics of entan- 
glement, remarkably similar, are characterized by a moder- 
ate growth. As a result, the evolution can be simulated, with 
reasonably modest computational resources, for times long 
enough that the perturbation, initially localized, can spread 
over large regions, involving e.g. of the order of one hundred 
sites. 

We notice that the TEBD algorithm has also been applied to 
other classes of evolution. In the case of a global quench that 
modifies the Hamiltonian everywhere in the chain, or when 
the initial state of the system does belong to the low energy 
sector of H, it has been reported that a large amount of en- 
tanglement is created lfl6l fl7ll and therefore the simulation 
becomes inefficient. We emphasize that this is not at all sur- 
prising, given the exponentially large dimension of the Hilbert 
space in the system size, which implies that a generic state of 
the system cannot be represented efficiently. Instead, what is 
truly remarkable (and useful for practical simulations) is that 
there are circumstances of interest where the state of the sys- 
tem can be efficiently encoded using an MPS. Apart from the 
case of a local perturbation of the ground state/local quench, 
there are not many other classes of evolution known that can 
be efficiently simulated for reasonably large times. One no- 
ticeable exception is the case of a global quench in a disor- 
dered system lfl7ll . which seems to be related to localization 
effects experienced by information in the presence of disorder 
frill . Another case, perhaps less surprising but still poten- 
tially very useful, is the evolution of observables and of ther- 
mal states after a global quench in integrable systems [19]. 

The rest of the paper is organized as follows. In Section II 
we introduce the measure of entanglement that is relevant in 
the present context and describe the one-dimensional models 
we have considered. In Section III we present the results of the 
simulations. Specifically, we describe two cases that illustrate 



the two types of scaling of entanglement observed in all the 
simulations we have performed. In Section IV we discuss the 
previous results, including simple toy models that reproduce 
the entanglement scaling reported in Section EH, and present 
some conclusions. 



H. ENTANGLEMENT MEASURE AND 
ONE-DIMENSIONAL MODELS 

Our goal is to characterize the growth of entanglement dur- 
ing the time evolution that takes place after the ground state 
l^cs) °f a ID system has been perturbed by some local opera- 
tor A, Eq. (Q~|i. In the present context, an appropriate measure 
of the entanglement contained in is provided by the 

number \ °f terms in its Schmidt decompositions (we refer to 
IH [H for details). More specifically, given a partition of the 
chain into two blocks containing the first r sites and the re- 
maining N — r sites respectively, the Schmidt decomposition 
reads 

i*w> = f>;i<i>r>i<^ +1:JV }, o) 

a=l 

where the rank x r > me Schmidt coefficients \ r a and the 
Schmidt bases | < l ) ^ r ) and |$q +1:JV ) are all time-dependent. 
Eq. (f3]) contains in principle a very large number of terms, 
with \ r ~ exp(TV). However, when the Schmidt coefficients 
decay fast with a, a good approximation to \^(t)} may 
be obtained by keeping only a relatively small number \l of 
terms. The truncation introduces a small error e l2(ill . given by 

e^l-J2(K) 2 - (4) 

The exact and approximate ranks \ and Xt of |\l/(t)) are then 
defined by maximizing over bipartitions, 

X = max x r \ Xe = max Xe- ( 5 ) 

r r 

A MPS can store the truncated state using 0(NxT) coef- 
ficients, while the cost of simulating a time step scales as 
O(Nxe)- Therefore Xc is indicative of the computational re- 
sources required in an approximate simulation with truncation 
error e. Although in an actual simulation there might be other 
sources of errors, such as those due to Trotter expansion (H, 
and errors accumulate in time, the simple entanglement mea- 
sure Xe turns out to be informative enough as to allow us to 
assess in practice whether a simulation is efficient. 

We have considered a large number of quantum models on 
a one dimensional lattice. Our hope is that by studying these 
models, we already observe all possible forms of scaling of 
entanglement, as measured by % e , so that we can draw con- 
clusions that are valid for generic ID systems — or at least 
for those ID systems that are usually of interest in condensed 
matter physics and quantum statistical mechanics. The spe- 
cific models we have considered are: 
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1. — The quantum Ising model with parallel and transverse 
magnetic fields 

II = E + E ( h *< + h *<) , (6) 

r r 

where h x and h y are the intensity of uniform magnetic fields 
in the x (parallel) and z (perpendicular) directions. For h x — 
0, a Jordan-Wigner transformation maps this quantum spin 
model (with spin |) into a model of free spinless fermions. 

2. — The quantum XY model with transverse magnetic field 

H " = E (^x +1 + ^r a X +1 + K<) , (7) 

where 7 is the anisotropy parameter and h r z the intensity of a 
(possibly site-dependent) transverse magnetic field. For 7 = 1 
we recover the quantum Ising model with transverse magnetic 
field, whereas 7 = corresponds to the quantum XX model, 
which is also used as a model of hard-core bosons. Again, a 
Jordan-Wigner transformation maps this quantum spin model 
(with spin |) into a model of free spinless fermions. 

3. — The quantum XXZ model with magnetic field in the z di- 
rection 

#*xz = E + « +1 + + K<) , (8) 

r 

where A is the anisotropy parameter and h r z the intensity of 
a (possibly site dependent magnetic field in the z direction. 
For A = we recover the XX model, for A = 1 a spin- 
i Heisenberg quantum antiferromagnet and for A = — 1 a 
model locally equivalent to a spin-i Heisenberg quantum fer- 
romagnet. Through a Jordan-Wigner transformation, H xxz be- 
comes a model of interacting spinless fermions. 

4. — The Bose-Hubbard model 

H BU = E [- J(b ts b s+1 + he.) - fm s + Un s (n s - 1)] , 

S 

(9) 

a model of interacting bosons where J is the hopping ampli- 
tude, [i is the chemical potential, U is the on-site repulsion; 
and 6^ s , b s and n s = b' s b s are the bosonic creation, annihila- 
tion and number operators at site s, respectively. 

For all the above models we have considered gapped and 
gapless cases. For the XY and XXZ models we have consid- 
ered, in addition to uniformed magnetic fields, disordered ver- 
sions by introducing inhomogeneous, random magnetic fields 



HI. RESULTS 

We have simulated time evolutions of the form (Q]) for each 
of the variants of the models ©-(O of the previous section, 
for chains with open boundary conditions and with a number 
TV ranging from 100 to 200. After computing the ground state 
|\T/ GS ) of a given Hamiltonian H through a simulation of imag- 
inary time evolution |3|. at / an operator A is applied to 



the middle of the chain. The system is then allowed to evolve 
according to the time evolution operator e~ lHt , that results in 
a propagation of the perturbation from the center of the chain 
toward its ends. We stop the simulations before the signal 
reaches the boundary of the chain. Chains of different lengths 
have been simulated in order to guarantee that our results are 
essentially independent of the system's size. Convergence of 
the results in this sense is only possible in those models where 
the correlation length is sufficiently smaller than the system's 
size. We notice that systems at or very close to a quantum crit- 
ical point fail to fulfill this condition. Also, the simulations 
have been performed with an MPS of rank x* much larger 
than the reported Xe> m order to guarantee that our analysis 
is independent of the computational resources devoted to the 
numerical simulation used to obtain them. Specifically, we 
report results for \e < 35 which have been obtained by simu- 
lating the evolution with an MPS with x* = 80, 100 and 150. 

The evolution in all simulations has been carried on with a 
fourth order Trotter expansion (labeled as Z\ in i22ll ) with a 
time step St = 0.02. The associated Trotter error is e Tloll oc 
St 5 x T/5t = 8 x 10~ 6 for T = 50 time units. The computa- 
tion times varied from few hours to three weeks in the worst 
case, in a PC with dual core processor at 2.2GHz with 4GB of 
RAM. In the Bose-Hubbard model, each site was truncated to 
four levels. 

We observe that, in spite of the rather diverse nature of the 
physics described by the models ©-(O, the scaling of its en- 
tanglement in time follows a very similar pattern. Next we 
describe in detail the results obtain for two specific systems. 
They have been chosen as representatives of the behavior of 
Xe in all the models under study. 

Example I: Saturation. — Our first example is concerned 
with a system of interacting bosons described by the Bose- 
Hubbard model H BU of Eq. (O in a chain with N = 100 sites. 
Fig. (QJ) shows the reorganization of the density of bosons 
when one extra particle is introduced in the middle of the 
chain, which is in the Mott insulating phase. Fig. (f2^) presents 
the changes in the spectrum of the squared Schmidt coeffi- 
cients (A£r 50 ) 2 during the evolution. We depict the Schmidt 
coefficients for r — 50 (middle of the chain) because this is 
the bipartition that appears most entangled (i.e., Xt = xl =5 ° 
for all e) at all times (except some minor corrections of one or 
two sites). We measure time in 1/U units, and choose U = 1. 

We notice that there are two marked regimes. First, for 
times up to t = 20 — 30, a number of Schmidt coefficients in- 
crease monotonically, with a growth that progressively slows 
down. Bipartite entanglement between the left and the right 
halves of the chain is being created (see an analogous behav- 
ior obtained analyticaly in Ref. lfl5tl ). Then, for larger times, 
and coinciding with the fact that the fronts of the density wave 
are far from the middle of the chain, the Schmidt coefficients 
become roughly stationary, indicating that the production of 
bipartite entanglement at the center of the chain has come to 
a halt. Notice that this is reflected in an initial growth, then 
saturation, of Xe> as depicted in Fig. (fJJ}). Elsewhere in the 
chain, say r = 30, no bipartite entanglement is generated un- 
til the front of the density wave arrives. At that point, the 
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Figure 1 : Propagation of a perturbation introduced in the middle of 
a chain of N = 100 sites. This evolution corresponds to a Bose- 
Hubbard model, Eq. in the Mott-insulating phase (hopping term 
J = 0.1, on-site repulsion U — 1 and chemical potential fj, = 0.5) 
with one boson per site, (a) Expectation value (n r ) for the number 
operator on each site, as a function of time, after an extra boson is 
introduced on site so = 50. The perturbation propagates throughout 
the chain, occupying a region that grows linearly in time until reach- 
ing the boundaries, (b) The entanglement entropy S r shows how 
every bipartition remains almost disentangled until the arrival of the 
wave front. When the front is gone far away, the entropy saturates. 

Schmidt coefficients start increasing to later become sta- 
tionary, in a pattern that mimics what occurred at the center of 
the chain, but with increasingly delayed and attenuated growth 
as we move away from the center. The entanglement entropy 
of a bipartition, 

S^-]T(A^ 2 log 2 (A^ 2 , (10) 

a 

offers a complementary, coarse-grained picture into the entan- 
glement growth in the system, that confirms the above obser- 
vations, see Fig. (QJ>). 

Example II: Moderate steady growth. — A slightly dif- 
ferent pattern of entanglement growth occurs in the antifer- 
romagnetic quantum Ising model with tilted magnetic field 
H Istag of Eq. (|6]l with N = 200 spins. Fig. (fJk) shows the 
propagation of a magnetization wave produced by flipping 
one spin in the center of the chain, r = 100, for the case 
h x = h z = 1. Fig. (0^) shows that, as in the first example, 
a number of Schmidt coefficients first grow substantially and 
then tend to saturate. An important difference, however, is 
revealed by studying Xe> Fig. &>)■ For e = 0.1, 0.01, the ap- 
proximate rank \ t again saturates, but for smaller e it keeps 
growing in time. This is due to the appearance of an increas- 
ingly large number of small Schmidt coefficients, that need 
to be consider collectively in order to decrease the truncation 
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Figure 2: Entanglement growth at the center of the chain for the 
Bose-Hubbard model, Eq. ([9), in the same scenario as Fig. (Q3. (a) 
Starting from a barely entangled state, with only a few non-vanishing 
Schmidt coefficients, several other coefficients start to grow, then sat- 
urate, (b) The approximate rank Xe first grows, then saturates at a 
number that increases as we decrease the truncation error e. This im- 
plies that an accurate approximation to | VP (t) ) can be stored in a MPS 
with a small rank Xe that after some time becomes roughly constant. 
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Figure 3: Propagation of a perturbation introduced in the middle of 
a chain with TV = 200 sites, corresponding to the Ising model with 
tilted magnetic field, Eq. (|6), with h x = h z = 1. Starting with a 
slightly entangled ground state, at t = one spin is flipped on site 
100 applying ctJ =10 °. The evolution of both (a) the expectation value 
(<Tz) for the magnetization on the z direction and (b) the entropy of a 
bipartition resemble their analogues in the Bose-Hubbard model, Fig. 
10. However, Fig. reveals that more entanglement is generated 
in this second system. 
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Figure 4: Entanglement growth at the center of the chain for the Ising 
model with tilted magnetic field, Eq. ® in the same scenario as Fig. 
I0. (a) As in Fig. l[2ji), starting from a barely entangled state, with 
only a few non-vanishing Schmidt coefficients, several other coeffi- 
cients start to grow, then tend to saturate, (b) But while rank Xe still 
saturates for e = 0.1, 0.01, it keeps growing (approximately linearly) 
for lower values e = 10 -3 , 10 -4 . [Notice that we have superposed 
the results of simulations performed with an MPS of increasing rank 
X* = 80, 100, 150. The results are generally very similar (except 
e.g. for a small variation in x e=10 -3)]. 



error e. Although it is difficult to characterize the scaling of 
Xe in time, we obtain results compatible with Xe ~ t p for a 
power p of the order of 1 . Correspondingly, a plot of the en- 
tropy, Fig. ((3J)), shows an initial rapid growth followed by a 
much slower, possibly logarithmic growth at longer times. 



IV. DISCUSSION 

Our simulations show a strikingly similar pattern of entan- 
glement generation in all the systems we have studied l2lll . In 
particular, we always find a sufficiently large e, e.g. e = 1CP 2 
in the examples of the previous section, for which Xe seems to 
saturate as a function of time. This means that an MPS with 
fixed rank Xe already suffices to approximate (to accuracy e) 
the state of the chain during the whole evolution. However, if 
we now try to improve this accuracy by lowering e, we find 
two possible behaviors. In some cases, as in example I, the 
rank x £ =io- 3 still saturates, although at a later time. In some 
other cases, as in example II, Xe=io- 3 keeps growing indefi- 
nitely (that is, as far as we could see with our computational 
resources), with a shape that resembles a small power of time 
Xe ~ t p . This moderate entanglement growth agrees with the 
analytical result obtained in Ref. Il2tl in a similar scenario. By 
means of conformal field theory they show that the entropy 
grows logarithmically with time, S(t) ~ log(f), when two 
initially decoupled quantum chains are joined together. This 
corresponds as well to a local quench of the ground state, and 
the subsequent evolution is expected to be analogous in terms 
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Figure 5: Toy model for the time evolutions of examples I and II 
of Sect. III. The chain is initially in a ground state with all spins in 
state |0). Either 1 or 2 spins are flipped into state jl) at the middle 
of the chain, and the evolution produces an entangled state |<6p) of 
Eqs. i ll It and d 1 2b (p = 1 or 2) that involves a number L of central 
spins that grows linearly in time. 



of entanglement creation through the whole chain. These mild 
behaviors are in sharp contrast with the exponential growth, 
Xe ~ e*, that one finds e.g., in a spin model where at time 
t = the value of the magnetic field throughout the whole 
chain is changed (global quench) or the spins in half of the 
chain are flipped id [M El- 

Although the scaling of Xe described in this paper corre- 
sponds to large times in a large systems, both the time and 
the system size are obviously finite. Our results show that 
at least during a valuably large period of time, the entangle- 
ment grows in a moderate way, enabling the efficient, accurate 
simulation of the evolution. Whether at even larger times the 
simulation remains efficient is hard to tell. This is not too rel- 
evant, however, in all those applications where the interesting 
phenomena occurs within moderate times from the moment 
the perturbation is introduced. 

In conclusion, we have studied the generation of entangle- 
ment in a particularly relevant class of time evolutions, namely 
those that follow from locally perturbing the ground state of a 
system, which are related to the computation of unequal-time 
two-point correlators (f2]i in one spatial dimension. Our results 
provide strong evidence that such evolutions can quite often 
be efficiently simulated with algorithms, such as the TEBD, 
based on updating a MPS |2|, [3l U, |a [3, l8[] . Our expectation is 
that the scaling of entanglement that we have observed for the 
specific models under consideration also applies to most one- 
dimensional models of interest in condensed matter physics 
and quantum statistical mechanics. It would be extremely in- 
teresting to be able to formally justify these results. This does 
not seem to be an easy task. 

We conclude by describing two simple toy model evolu- 
tions that reproduce the two observed patterns of entangle- 
ment scaling. Initially, a spin chain is in a ground state where 
all the spins are in state |0). In one case, the perturbation flips 
one single spin into state |1), and the time evolution produces 
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an entangled state 

\^i) = ^^\0i---U---0l) (ID 

involving the L central spins of the chain, where L grows lin- 
early in time, L = 2vt, with v is the speed of propagation of 
the perturbation, see Fig. (0. Here, |$f ) is a linear combina- 
tion of all the strings of L spins containing one single one. It 
can be seen that in any link of the chain, the Schmidt rank \ 
will go from 1 (before the expanding perturbation has reached 
that link) to just 2. Thus, this model reproduces the saturation 
observed in the example I. In the second case, the perturbation 
flips two spins into state |1) and the time evolution produces 
an entangled state 

= rrrr-T S E |0i---V--i 3 ----0l) (12) 



involving the L central spins of the chain, where as before L 
grows linearly in time, L = 2vt. Now |$2 ) is a linear com- 
bination of all the strings of L spins containing two ones. It 
can be seen that in any link of the chain, the Schmidt rank 
X grow linearly in time from the moment the expanding per- 
turbation reaches that link. Therefore this second evolution 
reproduces the moderate steady growth observed in the exam- 
ple II. Perhaps then (more sophisticated versions of) these toy 
models will inspire new simulation algorithms, based on an 
alternative ansatz that is even more efficient than an MPS for 
the considered class of time evolution. 
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